; docformat = 'rst'
;
;+
;
; :Purpose:
;   Calculate fractional array indices for use with interpolation routines
;
; :Inputs:
;   old: the "original" 1D array
;   new: the "new" 1D array
;
; :Outputs:
;   A 1D floating point array of fractional array indices
;
; :Example::
;   old=indgen(3)
;   new=[1.5, 3.3]
;   ii=mr_interp_indices(old, new)
;   print, ii
;     1.5     3.3
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, Mar 12, 2014
;
;-
function mr_interp_indices, old, new

  compile_opt idl2, hidden

  old=reform(old)
  new=reform(new)

  if n_elements(new) gt 1 then begin
    if sign(old[1] - old[0]) ne sign(new[1] - new[0]) then stop
  endif

  newsize=n_elements(new)
  oldsize=n_elements(old)

  index=dblarr(newsize)

  For NewI=0, NewSize-1 do begin
    wh=where(Old lt New[NewI])
    if wh[0] eq -1 then begin
      ;smaller than first element
      Index[newI]=0.d + (new[newI] - old[0])/(old[1] - old[0])
    endif else begin
      if n_elements(wh) eq oldsize then begin
        ;larger than last element
        Index[newI]=double(OldSize-1) + (new[newI] - old[OldSize-1])/(old[OldSize-1] - old[OldSize-2])
      endif else begin
        ;in between two grid cells
        index[newI]=double(max(wh)) + (new[newI] - old[max(wh)])/(old[max(wh)+1] - old[max(wh)])
      endelse
    endelse
  Endfor
  
  return, index

end